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Fundamental understanding of complex dynamics in many-particle systems on the atomistic level 
is of utmost importance. Often the systems of interest are of macroscopic size but can be partitioned 
into few important degrees of freedom which are treated most accurately and others which constitute 
a thermal bath. Particular attention in this respect attracts the linear generalized Langevin equation 
(GLE), which can be rigorously derived by means of a linear projection (LP) technique. Within 
this framework a complicated interaction with the bath can be reduced to a single memory kernel. 

This memory kernel in turn is parametrized for a particular system studied, usually by means of 
time-domain methods based on explicit molecular dynamics data. Here we discuss that this task is 
most naturally achieved in frequency domain and develop a Fourier-based parametrization method 
that outperforms its time-domain analogues. Very surprisingly, the widely used rigid bond method 
turns out to be inappropriate in general. Importantly, we show that the rigid bond approach leads 
to a systematic underestimation of relaxation times, unless the system under study consists of a 
harmonic bath bi-linearly coupled to the relevant degrees of freedom. 


INTRODUCTION 

Studying complex dynamics of many-particle systems 
has become one of the main goals in modern molecu¬ 
lar physics. The fundamental understanding of the un¬ 
derlying microscopical processes requires the interplay of 
elaborate experimental techniques and sophisticated the¬ 
oretical approaches. Experimentally, (non-)linear spec¬ 
troscopy revealed itself as a powerful tool for prob¬ 
ing the dynamics and for determining the characteristic 
timescales, such as dephasing/relaxation times and reac¬ 
tion rates to name but two. For interpreting the experi¬ 
mental spectra theoretical models are needed which can 
give insight into the atomistic dynamics. Often, a reduc¬ 
tion of the description to few variables is convenient in 
many cases since this can not only ease the interpretation, 
but enable the identification of key properties [T] . Such 
a reduced description can formally be obtained from the 
so-called system-bath partitioning, where only a small 
subset of degrees of freedom (DOFs), referred to as sys¬ 
tem, is considered as important for describing a physi¬ 
cal process under study. All the other DOFs, referred 
to as bath, are regarded as irrelevant in the sense that 
they might influence the time evolution of the system 
but do not explicitly enter any dynamical variable of in¬ 
terest. Practically, such a separation is often natural, for 
instance, when studying a reaction with a clearly defined 
reaction center or solute dynamics in a solvent environ¬ 
ment. Further, reduced equations of motion (EOMs) for 
the system DOFs can be derived in which the influence 
of the bath is limited to dissipation and fluctuations. 

The most simple formulation of this idea is provided by 
the Markovian Langevin equation, where dissipation and 
fluctuations take the form of static friction and stochastic 
white noise, respectively US!. Situations where memory 
effects become important are accounted for by the gen¬ 


eralized Langevin equation (GLE) [MZ] via a frequency- 
dependent friction and a stochastic force with a finite 
correlation time. This generalized equation has been em¬ 
ployed for instance, in the theory of vibrational relaxation 
for estimating characteristic relaxation times [SHII], re¬ 
action rates [12] and for thermostatting purposes [ISHISI; 
see e.g. Refs. msi] for review. 

The microscopic origin of the GLE can be rational¬ 
ized starting from different standpoints. First, it can 
be rigorously derived from the so-called Galdeira-Leggett 
(GL) model, where the environment is assumed to be a 
collection of independent harmonic oscillators bi-linearly 
coupled to the system jU [HI This model has 

been widely used in analyzing and interpreting (non- 
)linear spectroscopic experiments on systems in con¬ 
densed phase, termed multi-mode Brownian oscillator 
(MBO) model in this context [2TI - f26] . The second, more 
formal ansatz is to employ projection operator techniques 
in order to recast the system’s EOMs into linear or non¬ 
linear GLE forms UlillTlllT]. In the former, the resulting 
system potential is effectively harmonic, whereas in the 
latter the system potential is formed by a (non-linear) 
mean-field potential. In this approach noise and dissi¬ 
pation can be mathematically defined as (non-)linearly 
projected quantities. 

In any case, practical use of the GLE can only be made 
in connection with a stochastic model for the noise term 
being the main assumption of the formalism. The gen¬ 
eral advantage of the stochastic GLE is that dissipation 
and the statistical properties of the noise are entirely de¬ 
scribed by the so-called memory kernel being simply a 
function of time. If such a memory kernel can be obtained 
for a real system then the full quantum-mechanical treat¬ 
ment of the bath can be performed analytically, lead¬ 
ing to a quantum version of the GLE |28ll32j . Alterna¬ 
tively, a density matrix theory either via the Eeynman- 
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Vernon influence functional approach |331134] or hierar¬ 
chy type EOMs [35l [36] can be employed. Further, nu¬ 
merical methods that solve the Schrddinger equation in 
many dimensions m or various quantum-classical hy¬ 
brid schemes [38] can be used. Finally, the machinery for 
a purely classical treatment of the GLE is provided by 
the method of colored noise thermostats mils] or simi¬ 
lar techniques [39l |40] . Therefore the GLE formalism has 
become a popular tool for assigning system properties in 
macroscopic environments. 

However, establishing a connection between a real 
molecular system and a GLE might not be straightfor¬ 
ward. In a recent study we have shown that in condensed 
phase the mapping between the two can be established 
only for the effectively harmonic GLE derived by means 
of linear projection (LP) operator techniques [H]. All 
other possible mappings onto the GLE where the system 
force is kept anharmonic turned out inapplicable due to 
the so-called invertibility problem. Hence, in this work 
we limit ourselves exclusively to parametrizing the linear 
GLE. 

Common approaches to calculating memory kernels 
involve classical molecular dynamics (MD) simulations 
where the environment is explicitly taken into account. 
A very popular scheme is to obtain the memory kernel 
from the time-correlation function (TCF) of the forces 
exerted on a frozen system coordinate, referred to as 
the rigid bond approach 0 isi in HU- However, Berne 
et al. [33] have shown that this ansatz is only correct 
when the system frequency is much larger than the bath 
ones. If such a frequency separation is not given, one can 
determine the memory kernel from a Volterra integro- 
differential equation for the momentum-autocorrelation 
function (MAE). Practically, the memory kernel can be 
computed from explicit MD MAFs involving discretiza¬ 
tion schemes in time domain [44ll47] . An example of this 
ansatz is the method introduced by Berne and Harp who 
have employed polynomial interpolations of the MAE in 
order to calculate necessary derivatives [iniiiiiiiiHj- 

Another idea is based on Laplace domain tech¬ 
niques [39l SSI |49| . Here, the Volterra integro-differential 
equation is transformed into Laplace domain, resulting 
in an algebraic expression for the Laplace-transformed 
memory kernel. However, the practical application of 
this method has not been discussed in detail. In contrast, 
this technique has been criticized due to the significant 
difficulties of the numerical Laplace back transform [47] . 

In this paper we demonstrate that it is more convenient 
to parametrize the GLE in frequency domain, where a 
memory kernel turns into a spectral density. We present 
a new Fourier-based method, which provides a direct and 
robust way for calculating spectral densities of realistic 
solute dynamics in liquid solvents and avoids the afore¬ 
mentioned numerical problems. It is shown that exist¬ 
ing time-domain techniques suffer from either numerical 
or conceptual dehciencies. Very surprisingly, the well- 


established rigid bond approach turns out to be inappli¬ 
cable, even when the frequency separation is provided. 

The paper is organized as follows. After this intro¬ 
duction, the GLE formalism is described in detail. Then 
the rigid bond approach and the method by Berne and 
Harp, [44] which are widely used for calculating the mem¬ 
ory kernel from explicit MD simulation data are reviewed. 
Further, we present in detail a new efficient scheme that 
enables calculating the memory kernel in Fourier space. 
Finally, we demonstrate that our procedure gives reason¬ 
able spectral densities for realistic solute dynamics on the 
example of two hydrogen bonded systems: HOD in H 2 O, 
and the ionic liquid [G 2 mim][NTf 2 ]. The model system 
introduced by Berne et al. [43], which consists of a di¬ 
atomic in an atomic gas, termed A 2 in A, is considered for 
cross-checking. A comparison of our procedure against 
the rigid system approach and the method of Berne and 
Harp is provided, and the failure of the rigid bond ap¬ 
proach is analyzed in detail followed by conclusions and 
outlook. 


GENERALIZED LANGEVIN EQUATIONS 

In the following we consider a classical one-dimensional 
system with a mass m, coordinate x and conjugate mo¬ 
mentum p, which undergoes Brownian motion in a clas¬ 
sical bath described by coordinates {Qi} and momenta 
{Pi}. The total open system is without any loss of gen¬ 
erality assumed to be partitioned into the system, Vg(a;), 
and the bath, V^{{Qi}), coupled via a system-bath in¬ 
teraction, Vg_B(a:, {Qi}). 


The linear GLE 

A mathematically rigorous approach for deriving re¬ 
duced EOMs is to employ LP operators, which project 
a dynamical variable A{x,p', {Qi, Pi}) depending on the 
full set of system and bath coordinates onto the linear 
subspace spanned by x and p. mum The resulting lin¬ 
ear GLE reads 


P{t) 

x{t) 


kT 


t 

At) - J ^Lpit - 'r)p(r)dr R{t) 


Pjt) 

m 


( 0 . 1 ) 


where T stands for the temperature, k is the Boltzmann 
constant and (...) denotes canonical averaging. The dis¬ 
sipative force is thus given by a convolution of the mo¬ 
mentum with a memory kernel, ^lp(O) which is a real 
function decaying on a finite timescale. The noise term, 
R{t), formally contains an explicit dependence on all 
bath DOFs and is related to the memory kernel via the 
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fluctuation-dissipation theorem (FDT) 

(i?(0) exp[iLLpt]i?(0)) 
- 


spectroscopy pTfl26) . This model assumes that the bath 
consists of independent harmonic oscillators bi-linearly 
coupled to the system. The total potential energy for 
the model reads 


where ipp is the Liouville operator that describes the 
linearly projected time evolution; note that it does not 
correspond to the real Hamiltonian flow. In practical 
applications one sacrifices the deterministic time evolu¬ 
tion of the noise and mimics this term by a stochastic 
zero-centered Gaussian process R(t) keeping the FDT as 
its main statistical property [T51 inni HO]- The GLE in 
Eq. (0.11 then becomes a non-Markovian stochastic dif¬ 
ferential equation, where memory effects in the noise are 
incorporated via the FDT, Eq. (0.2). The Markovian 
limit Clp(^) ^Q5{t) can be obtained if the decay of 
the memory kernel is much faster than a characteristic 
timescale, e.g. a vibrational period, of the system. The 
great advantage of the stochastic GLE is that the im¬ 
plicit characterization of the bath is simply given by the 
memory kernel Clp( 0- Tor rm easier interpretation it is 
convenient to consider ^lp(^) in frequency domain, where 
it is referred to as the spectral density 


V{x,{Qi}) = Vs{x) -'^QiQ^x (0.5) 

i i 

with the bath frequencies Wj, bath masses set to unity and 
the coupling strengths g,. [laiTiin] Often the square 
on the right hand side of Eq. ( |0.5[ ) is completed causing 
thereby a harmonic correction to the system potential 
Vs{x) 

Vsix) = Ts(a:) - i ^ %x^ . (0.6) 

I * 

The corresponding GLE implied by the CL model can 
be derived both in classical Hi] and quantum nnin] 
domains. Limiting ourselves to the classical description, 
the derivation can be straightforwardly performed with¬ 
out any further approximations by integrating the EOMs 
for the bath, yielding 


Ilp(w) = J exp[-iuji\^i^p{t)dt . (0.3) 

0 


p(0 = -^-/^Chit - x)p{T)dT + R{t) ; (0.7) 

0 


Here and in the following hat denotes the half-sided 
Fourier transform. It should be stressed at this point 


that the term which is linear in x in Eq. (0.1), usually 


referred to as the mean intramolecular force of the sys¬ 
tem, [M El is a consequence of the LP. This implies 
that one can interpret the projected system part as an 
effective harmonic oscillator with the frequency 


-2 _ 
U! = 


kT 

m{x^) 


Within the CL model the real part of the spectral density. 


Eq. (0.3), can be interpreted as the bath modes’ density 


of states weighted with the strength of the coupling to 
the system 


3fi4L(w) = E - ^i) 


( 0 . 8 ) 


(0.4) By comparing Eq. (0.8) and Eq. (0.6) one finds that the 


frequency renormalization term is given by 


even though the original system potential can be arbitrar¬ 
ily anharmonic. The anharmonicity is formally projected 
onto the bath and, hence, incorporated into the noise 
term and the memory kernel. Thus, using this model 
to disentangle, e.g. anharmonicity and mode coupling, 
by means of non-linear spectroscopies is doomed to fail 
from the outset. Still, if one is interested only in linear 
properties, like transport coefficients or linear absorption 
spectra, the GLE gives the correct description im and 
provides a very simple and intuitive model. In the gen¬ 
eral case the linear GLE does not reflect all dynamical 
features of the real system and, thus, one partly looses 
the true atomistic picture when using it. 


Ts(a:) = Vs{x)- — / 3?|cL(w)dw 


(0.9) 


This frequency renormalization term thus induces a red- 
shift of the system frequency due to interactions with 


the harmonic oscillator bath. Comparing Eq. (O.I) to 


Eq. ( |0.7| ) immediately suggests that the only possibility 
for them to coincide is to set the model system potential 
V^{x)=kTx^/{2{x^)). 


PARAMETRIZING A SPECTRAL DENSITY 


The Caldeira-Leggett Model 

Another way to derive a GLE is to employ the CL 
model that has enjoyed popularity in condensed phase 


Formulation of the problem 

In order to utilize the framework of a stochastic GLE 
one has to find a memory kernel that correctly mimics the 
features of the bath. We start the consideration with the 
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linear GLE, Eq. (0.1). Multiplying it with p(0), taking 
the canonical ensemble average and noting that the corre¬ 
lation function with the noise term vanishes by construc¬ 
tion, yields a Volterra integro-differential equation for the 
normalized MAE defined as C^pit) = (p(0)p(t))/(p^) 


t 

Cppit) =-mui^Cp^it)- J ^Lp{t - T)Cpp{T)dT , ( 0 . 10 ) 
0 

where the function Cpx{t) = {p{0)x{t))/{p^) denotes the 
momentum-position cross correlation. The aforemen¬ 
tioned equation can be simplified by using the relation 


Cpxit) = CppMdT , (0-11) 

0 

which yields 

t 

Cpp{t) = - J Ki^pit - r)C7pp(r)dr , (0.12) 

0 

with the shifted memory kernel 

^Lp(0 — Clp(^) + • (0.13) 


The basic idea is to invert this equation to obtain the 
memory kernel Kpp{t) and, hence, the spectral density 
Kpp{uj). As an input, the MAE calculated from an ex¬ 
plicit MD simulation is needed. 

Performing the aforementioned procedure with 
Eq. (0.7) as a starting point, yields another Volterra 
equation 


t 

Cpp{t) = Cppit) - J ^Chit - T-)Cpp{T)dT , (0.14) 

0 

where Cpp{t) is momentum-system force correlation 
function. Thus these two functions, Cpp{t) and Cpp(t), 
can be used to parametrize ^cl(^)- K should be stressed 
that a self-consistent procedure can be established only 
if the system is of the CL form. In general, there exist 
infinitely many pairs of the TCFs that correspond to the 
given spectral density, which is in the heart of the invert- 
ibility problem. m Still, the obtained spectral density 
is unique and can be used for comparison, see Sec. P 


is the derivative of Eq. (0.12) 


differential Eq. (0.12) itself. 


rather than the integro- 
Following Ref. [44], the 
First, the correspond- 
formulated by making a 


method consists of three steps, 
ing finite difference scheme is 
variable substitution r —>■ r — t in Eq. (0.12) in order 
to shift the t-dependence from 
differentiation with respect to t 

Kpp{mAt) = —Cpp{mAt) 


the kernel, followed by a 
resulting in 


-At^WjKpp{jAt)Cpp {{m-j)At) . (0.15) 

Here the time integration dTKpp{T)Cpp{t — t) is ap¬ 
proximated with the help of the Gregory formula employ¬ 
ing the weights Wj and MD timestep At. Note that the 
Gregory formula is advantageous to the Simpson rule, 
since it does not rely on the parity of the number of in¬ 
tegrated points [52] . 

Second, the MAE is interpolated by a 2nth order poly¬ 
nomial C'fc(t) = vicinity of fc-th time 

step such that Cfc((fc -|- j)At) = Cpp{[k j)At), for 
J = —n, —n + 1, ■ ■ ■ ,n — l,n. The derivatives that en¬ 


ter Eq. (0.15) are then approximated by the derivatives 


of the polynomial Cf^ at these time points. 

Third, one determines the initial values of Kpp at 
first four time instances. The first value ArLp(O) = 
(p^(0)/y(0)) is directly calculated from averaged MD 
data. Subsequently, the remaining starting val¬ 
ues Kpp{At),Kpp{2At),Kpp{3At) are accurately de¬ 
termined by applying the so-called Day’s method to 
Eq. dOT^ . [Sg 


Rigid Bond Approach 

Another time domain technique which has been widely 
used in the literature, e.g. in studies of vibrational relax¬ 
ation 01111113], is based on approximating Clp( 0 by 
^rb(^) calculated as 

&B(t) = ^ (0.10) 

Here, Lpg is the Liouvillian corresponding to the system 
with the frozen bond and the noise is calculated via ex¬ 
plicit MD simulations. If a system coordinate x is chosen 
as a single bondlength, as in the practical applications 
considered here, then the noise R{t) is obtained as 


Berne and Harp Method 


The method presented in this section approaches the 
iterative solution for Kpp{t) from Eq. (0.12) in the 
time domain has been introduced by Berne and Harp in 
1970 [H]; due to stability reasons, the starting point here 


R{t) = 


Flit) 

Till 


Ml 

m2 


(0.17) 


where Fi{t), F 2 {t) are the forces exerted on the two 
bonded atoms by their surrounding, , m 2 are their 
masses and ^ is the reduced mass of the atom pair. The 
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vector ni 2 is the unit bond vector pointing from atom 1 
to atom 2. 

In general, fixing the system’s coordinate has an influ¬ 
ence on the energy flow between system and bath thereby 
affecting the memory kernel. The consequences have 
been extensively discussed by Berne et. al [23] and are 
briefly summarized below. In general, there are two ap¬ 
proximations behind the rigid bond approach, which can 
be formulated in terms of the memory kernels 'fLp(^)i 
defined in Eqs. (0.2 0.161 and the time corre¬ 
lation function due to the Hamiltonian flow of the real 
system 


^Rs(^) — 


(j?(0) exp (iLRs^)i?(0)) 
mkT 


(0.18) 


As it was shown by the authors, a simple relation between 
s-iid ^Lp(^) can be obtained in Laplace domain 


Irs(s) = - 


Ilp('S) 


1 -I- s/(w^ -I- s"‘)^Lp(s) 


(0.19) 


where tilde denotes the Laplace transform throughout 
the manuscript. Taking the limit a) —>■ oo in Eq. (0.19) 
naturally results in lim^^^ Irs(s) = Ilp(s)- As 

a matter of fact, the limit also directly corresponds to the 
rigid bond dynamics, fRB(s) = lim^_,.oo |rs(s). This can 
be most easily understood in the time domain, namely, if 
the memory kernel varies slowly on the timescale of the 
fast oscillations of the system coordinate, then the dissi¬ 
pative term in the GLE would vanish, as it is required by 
the rigid bond approach. Combining the two expressions 
obtained in the limit w —)■ oo yields the aforementioned 
two approximations ^lp(^) ~ Crs(^) ~ Crb(^)) see Sec. | 
for further discussion. 

Importantly, if the system studied would indeed be of 
the CL form, then the match ^cl(0 = ^rb( 0 becomes 
exact and thus one can exclude Crs( 0 from considera¬ 
tion. The memory kernels Ccl( 0 ^cid ^lp(^) coincide ex¬ 
clusively for harmonic system potential. Vs (a:), see Sec.[l 


Motivation for a Fourier Method 


The aforementioned methods are exclusively time do¬ 
main techniques. Another possibility that has been 
suggested in the literature is based on a transform of 
Eq. (0.12) into Laplace domain. |31||231I21] This ansatz 


results in an algebraic equation 


sCppis) - Cppit = 0) = -iLLp(s)Cp„(s) . (0.20) 


Note that a differentiation in time domain results in a 
multiplication by s in Laplace domain and the convolu¬ 


tion in Eq. (0.12) thereby turns into a simple product. 


Despite its apparent simplicity, the detailed implementa¬ 
tion of the method has not been carried out to our best 
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FIG. 1: Fourier transform (left panel) and Laplace transform 
(right panel) of a memory kernel consisting of a superposition 
of two exponentially damped cosine functions. Their individ¬ 
ual Laplace/Fourier counterparts are plotted in black dotted 
lines. 


knowledge, owing to numerical instabilities of a Laplace 
back transform. m- Further we show that both time and 
Laplace domain methods are not natural for the present 
purpose. 

First, numerical algorithms for GLE simulations of¬ 
ten require to fit the memory kernel to a special class 
of analytic functions, typically exponential and/or expo¬ 
nentially damped cosine functions [aiini [3i unuMj- 
Practically, a fit to oscillatory functions in time domain 
can be very difficult since the memory kernels of realis¬ 
tic solute systems usually constitute a mixture of terms 
with distinct frequencies. Importantly the same problem 
occurs in Laplace domain, see panel b in Fig. thereby 
making a fit of the memory kernel directly in Laplace do¬ 
main not recommendable as well. In contrast, the signals 
can be well separated in frequency domain, see panel a 
therein, thus simplifying the fit remarkably. This sug¬ 
gests that a successful method might be formulated in 
the frequency domain. 

Another benefit of working in frequency domain is the 
possibility to limit the description to the region of the 
spectral density that is resonant with the system fre¬ 
quency only. This is illustrated in Fig. by comparing 
the linear absorption spectra of a harmonic oscillator in 
a bath described by a spectral density with and without 
off-resonant contributions. The two spectra are almost 
identical and the influence of the off-resonant peaks in the 
spectral density on spectra is negligibly small, see insets 
therein, although the off-resonant peak intensity is about 
five times larger than the intensity in the resonant region. 
This is a manifestation of the well-known fact that the 
energy can be exchanged efficiently only when the bath 
is resonant with the system. Thus the consideration can 
be limited to a narrow frequency interval around the sys¬ 
tem frequency. In practice, this is a great simplification, 
since the spectral density might have an extremely elab¬ 
orate pattern, if the bath consists of fairly complicated 
molecules, see, e.g. panel c) in Fig. In time domain 














6 



m (arb. u.) 


FIG. 2: Two spectral densities, one having only the resonant 
contribution (green) and the other having additionally two off- 
resonant Lorentzian contributions (red) are depicted in the 
upper panel. The corresponding spectra of the harmonic os¬ 
cillator with frequency oj = 0.4 coupled to the spectral den¬ 
sities are shown in lower panel. Insets zoom on the spectral 
features stemming from the off-resonant contributions in the 
second spectral density. 


it would be impossible to filter out the corresponding 
signal. Importantly it would be equally complicated to 
extract the frequency window in Laplace domain due to 
the unsuitable shape of the Laplace-transformed damped 
oscillations, see panel b in Fig. 

To this end the problem of parametrizing GLE simula¬ 
tions seems to naturally pose itself in frequency domain, 
thereby operating with the spectral density rather than 
with the memory kernel. The transition from Laplace 
to frequency domain can be easily achieved by setting 
the Laplace variable imaginary, i.e. s = iui with a real 


frequency oj. Equation (0.20) then becomes 


iuCM - C„Jt = 0) = . (0.21) 


Solving Eq. (0.21) for the spectral density yields 

1 




( 0 . 22 ) 


-'ppV 


using that Cpp(O) = 1 by normalization; note that the 
convolution theorem still applies to the half-sided Fourier 
transform. 


Equation (0.22) is the starting point of the proposed 


Fourier-based GLE parametrization scheme, later re¬ 
ferred to as the Fourier method, that is presented in detail 
in the next section. 


Implementation of a Fourier Method 


In this section we elaborate on setting up a GLE simu¬ 


lation according to Eq. (0.1) using the previously derived 


result in Eq. (0.22) as a starting point. Two quantities 


need to be calculat ed: the effective harmonic frequency, 
a), defined in Eq. (0.4) and a spectral density ^Lp(a;). 


Since the two memory kernels ^lp( 0 iLLp(0 differ 
only by the constant the corresponding half-sided 
Fourier transforms are connected via 

-2 

iTLp(w) = |lp('^) +— f— . (0.23) 

UJ 


One notices that the real parts of ^lp('^) Ki^p{u}) 
coincide for a; > 0. Since the memory kernel ^lp(^) is a 
real function, the knowledge of the real part of 'Clp('^) 
and, thus, of K-pp{uj) is sufficient, because the imaginary 
part is provided by the Kramer-Kronig relations. Finally, 
fitting the hyperbola in the imaginary part of Kpp{oj) to 
the fit function h^co) = —d/uj provides an easy way to 
determine the effective harmonic frequency as d) = '/d. 

In order to use the obtained spectral density in GLE 
simulations, its corresponding memory kernel ■Clp( 0 
must be given as a superposition of damped cosine func¬ 
tions, f{t) = 2a^e~*'* cos(ct), ya,b,c S R, 6 > 0. The 
real part of its Fourier transform reads as a superposi¬ 
tion of 


3?/(a;) = a^b- 


1 

.6^ + (c — w)^ 


1 

b^ + {c + w)^. 


(0.24) 


which can in turn be used to fit the spectral density, 
^Lp(w). For fitting the hyperbola in the imaginary part, 
S5iTLp(a;), it is recommendable to subtract from it the 
imaginary parts of the fit functions given by 


3/(w) = 


C + UJ 


c — UJ 


+ {c + uj)^ + {c — uj)^ 


. (0.25) 


The performance of the method is illustrated on a sim¬ 
ple test system: a harmonic oscillator (w = 0.4) in a bath 
described by a spectral density comprised of a single func¬ 


tion as given in Eq. (0.24) with a = 0.03, b — 0.03 and 


c = 0.4. In order to test the procedure, the aforemen¬ 
tioned spectral density was used in GLE simulations to 
produce the MAF, that in turn was used to parametrize 
the spectral density as it was described above. The re¬ 
sulting spectral density was compared against the input 
one. 

It turns out that the successful use of the procedure 
presented above rests upon two important numerical is¬ 
sues. First, the Fourier transform Cpp{uj) needs to be 
calculated very accurately. In particular, performing the 
Fourier transform by means of the standard FFTW3 
library (based on a single-sided sum rule), led to 
diverging spectral densities even for moderate timestep 
sizes (0.1 — 0.3), that are typical for molecular systems 
with hydrogen atoms, if reasonable units (fs) are em¬ 
ployed, see Fig. Here, a more accurate quadrature, 
such as the 3/8 Simpson integration scheme, yielded ex¬ 
cellent results already for a time step of 0.3. 
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FIG. 3: Real part of spectral densities obtained according 
to Eq. (0.22 I for the test system, see text. The exact spectral 
density is shown in black. Numerical results obtained using 
the FFTW3 library for three time steps 0.1, 0.2, 0.3 are shown 
in red, green and yellow, respectively. The result employing 
the Fourier transform based on the Simpson rule for a time 
step of 0.3 is shown as blue dots. 



FIG. 4: Real part of spectral densities obtained for the test 
system, see text, upon smoothing the MAF. The exact spec¬ 
tral density is shown in black. Results for different window 
widths T = 50,100, 500, 5000 are shown as red, green, blue 
dots and yellow curves, respectively. 


The second numerical issue concerns the usually in¬ 
sufficiently converged tail of the MAF, that causes noisy 
results upon the Fourier transform. In order to reduce 
the noise level we suggest Gaussian filtering, that is to 
multiply the MAF by a Gaussian window function 


G{t) = exp 



(0.26) 


In frequency domain this corresponds to a convolution 
with a Gaussian function of the width Aw = 1/T thereby 
suppressing the noise. The parameter T should be chosen 
as a compromise between noise reduction and smoothing 


errors, see Fig. As a rule of thumb, T can be set equal 
to the correlation time in the system, which has been 
T = 500 for our test system. 

The developed parametrization method can be sum¬ 
marized by the following steps 


Calculate the MAF from explicit MD simulations 
with sufficient convergence (system dependent) 
Perform a reasonable Gaussian filtering to reduce 
the noise level upon the Fourier transform, taking 
the correlation time as a starting guess for the win¬ 
dow width 

Transform MAF into frequency domain using a suf- 
hciently accurate integrator, e.g. Simpson 3/8 rule 
Calculate the spectral density Ki^p{uj) using 


Eq. (0.22) 


Estimate w and fit the real part of Kpp{uj) in 
its vicinity to superpositions of functions given in 


Eq. (0.24) 


• Subtract the imaginary counterparts, Eq. (0.25) of 
the fit functions from 3ArLp(w) 

• Obtain the effective harmonic frequency w via a 
hyperbolic fit from the imaginary part of Kpp{uj) 

As a word of caution it should be stressed that the pre¬ 
sented scheme does not give a correct estimate of a spec¬ 
tral density at low frequencies due to the diverge nce o f 
Kpp{uj) at w = 0 caused by the J-function, see Eq. (0.23). 
This leaves description of pure dephasing times outside 
reach. 


MODELS AND COMPUTATIONAL DETAILS 

We have investigated two hydrogen-bonded systems - 
HOD in H 2 O, and the ionic liquid [C 2 mim][NTf 2 ] at room 
temperature, as well as an A 2 in A model system em¬ 
ployed by Berne et al. [43]. These examples have been 
chosen as water is perhaps the most important solvent 
and HOD features a spectrum with three distinct peaks, 
which makes it highly suitable for methodological investi¬ 
gations. [371 [55H55] The ionic liquid represents a system, 
where, in addition to moderate H-bonding, there exists 
a strong Coulomb interaction between the ion pairs [60] . 
Einally, investigating the A 2 in A model system allows 
us to compare our results with that of Berne et al, see 
Ref. [13] for the parameters. 

The aqueous systems have been comprised of 466 
molecules in a periodic box with the length of 2.4 nm 
interacting according to the force field adopted from 
Ref. [3T]. The harmonic HOD simulations used in 
Sec. [j for comparison, have been performed with ex¬ 
actly the same setup as the anharmonic ones, with the 
OH stretching potential in the HOD molecule being 
harmonic with the frequency obtained from expanding 
the respective Morse potential up to the second order 




























{^iuP' = 5081 kJ-mol~^A~^). The ionic liquid simula¬ 
tions have been carried out in a periodic box with the 
length of 4.5 nm containing 216 ionic pairs and the force 
field described in Ref. [52]. All explicit molecular dy¬ 
namics (MD) simulations have been performed with the 
GROMACS program package (Version 4.6.5) [53]. The 
spectra have been computed for the OH stretch in the 
HOD molecule and for the C(2)—H stretch of the imida- 
zolium ring |62j in the ionic liquid. In the latter case 
the potential for the C-H stretching motion has been 
re-parametrized to a Morse potential using DFT-B3LYP 
calculations [51]. Note that in the spirit of the system- 
bath treatment the 0-H bondlengths have been taken as 
the respective system coordinates, x, whereas all other 
degrees of freedom constituted the bath. We employed 
the “standard protocol” for calculating IR spectra, that 
is, a set of NVE trajectories, each 6ps long (time step 
0.1 fs), has been started from uncorrelated initial condi¬ 
tions sampled from an NVT ensemble. The dipole auto¬ 
correlation functions for the system coordinate have been 
Fourier-transformed to yield the spectra |65j . In order to 
achieve convergence, 1000 trajectories for the considered 
stretching motion have been employed. 

For GLE simulations we adopted the method of Gol- 
ored Noise thermostats [Til IT51151] . htting the spectral 
density to a superposition of 8-12 functions as given in 


Eq. (0.241. In the GLE simulations all the numerical pa¬ 


rameters for the time step, length, number of trajectories, 
etc. have been the same as in the explicit MD simulations. 
Bonds in the rigid bond approach have been hxed using 
the Settle algorithm as implemented in GROMACS. 




CL 
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FIG. 5: Real (left panels) and imaginary (right panels) parts 
of spectral densities (red curves) obtained according to the 
Fourier method are shown for a) A 2 in A, b) HOD in bulk 
water, and c) the ionic liquid. Vertical dashed lines denote 
the effective harmonic frequency uj. The fits to Lorentzian 
functions, Eq. (0.24|, performed in the resonant region and 


the hyperbolic fits to the imaginary parts are depicted with 
black crosses. The inset in panel c) zooms into the resonant 
region. Note the different scales for imaginary and real parts. 


NUMERICAL TESTS 

In this section we first demonstrate that the proposed 
Fourier method provides accurate spectral densities that 
can describe realistic solute dynamics. Then the perfor¬ 
mance of the time domain techniques is discussed, using 
the spectral densities obtained via the Fourier method as 
a reference. Finally, an exceptional case for which the 
rigid bond approach works, the A 2 in A model system, 
is discussed in detail. 


The Fourier method 

In order to test the Fourier method, first the explicit 
MD simulations are performed and vibrational spectra 
and MAFs are calculated. Then the spectral densities are 
parametrized from these MAFs according to the Fourier 
method. These spectral densities are used as input for 
implicit GLE simulations, which in turn yield vibra¬ 
tional spectra that are compared against their explicit 
MD counterparts. The coincidence of the GLE and MD 
spectra manifests the success of the Eourier method. 


In Fig. [^the spectral densities, ArLp(w)i calculated ac¬ 
cording to Eq. (0.221, are shown for the three systems 
studied. The real parts of ArLp(a;) are depicted in the 
left panels therein. For the A 2 in A model system, panel 
a), the spectral density nicely reproduces that obtained 
by Berne et al. [43] The spectral densities for the real 
systems considered, panels b) and c), possess peaked 
contributions at various frequencies stemming from the 
coupling to a plethora of vibrational modes in the bath. 
The spectral density is most structured in the ionic liq¬ 
uid case, panel c), due to the complexity of the molecules 
involved. However, as it was explained above, the only 
important region is in the vicinity of the system frequency 
w, denoted by a dashed vertical line therein. The fits to 
the real part of spectral density in the important region, 
black crosses in Fig. [^ illustrate the excellent fit quality. 
Finally, one sees that the hyperbola in the imaginary part 
(right panels therein), Eq. (0.231, is clearly pronounced 
and can be thus very well fitted. 


Having established the sufficient fit quality for the 
spectral densities, the GLE simulations are performed 
and the resulting vibrational spectra of the three in¬ 
vestigated systems (red stars) are compared with that 
obtained from explicit MD simulations (black curves) 
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FIG. 6: Explicit MD (black curve) and GLE spectra (red 
stars) shown for a) A 2 in A, b) HOD in bulk water and c) the 
ionic liquid. 


in Fig. It becomes apparent that the GLE simula¬ 
tions excellently reproduce the explicit MD spectra for 
all systems studied thereby confirming the success of the 
Fourier method as a parametrization scheme. 


Comparison to Time Domain Techniques 




500 1500 2500 3500 

CO (cm'^) 


Having obtained the numerical evidence that the spec¬ 
tral densities calculated from the Fourier method (red 
lines in Fig. are reliable, we compare the results ob¬ 
tained via the rigid bond (blue) and Berne and Harp 
methods (green) against them. One sees that the real 
parts of the spectral densities due to the Berne and Harp 
approach reproduce the reference results reasonably well. 
However, we encountered several numerical issues when 
using the method for the realistic systems studied. First, 
the polynomial interpolation turned out to be very sen¬ 
sitive to the degree 2n of the polynomial. In particu¬ 
lar, using n = I was insufficient, n > 2 led to notice¬ 
able spikes causing divergences in the kernel and only 
n = 2 yielded satisfactory results. Importantly, small 
inaccuracies in the polynomial interpolation turn out to 
accumulate strongly due to repeated integration accord¬ 
ing to Eq. (0.15). Second, even for n = 2 the kernel 
diverged at long times and one had to cut it manually 
at the plateau, which could not be determined uniquely 
(data not shown). These dehciencies together with the 
general problems of time domain methods discussed in 
Sec. [] make the use of Berne and Harp method not con¬ 
venient. 

Considering the results of the rigid bond approach 
(blue lines in Fig. [^, one sees that it is successful only 
for the A 2 in A system, panel a) therein. The curves 
corresponding to realistic systems deviate from the ref¬ 
erence results significantly. The general trend is that 
the larger is the frequency, the more the real part of the 
spectral density is downscaled. For HOD in water, panel 
b), noticeable deviations start from « 1000 cm”imply- 


FIG. 7: Real parts of the spectral densities obtained via the 
Fourier (red), rigid bond (blue) and Berne and Harp (green) 
methods are shown for a) A 2 in A b) HOD in bulk water and 
c) the ionic liquid. Dashed vertical lines mark lu position. 
Insets zoom on the resonant regions. 


ing wrong description in the region where the bending 
and OH stretching modes reside. Remarkably, in the 
case of the ionic liquid, panel c), the frequency region 
below 1500 cm” ^ is well reproduced by the rigid system 
approach. In contrast, contributions to the important 
resonant region are completely absent in both cases, see 
insets in Fig. [^for zoom on. Since the spectral densities 
obtained from the Fourier method are proven to be trust¬ 
worthy, we thus conclude that the rigid system approach 
breaks down for the realistic cases considered here. Ac¬ 
cording to the Laplace domain based analysis by Berne at 
al. [33], see also Sec.[j this breakdown is no surprise since 
the rigid approach is valid only if the system frequency 
is high compared to that of bath modes. However, as it 
was discussed in Sec. || it is more natural to consider the 
memory kernel in frequency domain. Here, the formula 
equivalent to Eq. (0.19|) reads 


?Rs(^) — 


Clp(w) 


1 -I- ioj /(w - a; )^Lp(a;) 


(0.27) 


From this equation it becomes clear that a problem 
emerges for uj ^ u), since the denominator diverges un¬ 
der the assumption that is finite in the vicinity 

of w; note this assumption is not fulhlled the irrelevant 
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case of isolated unperturbed systems. This divergence 
thus leads to vanishing ^rs(^) fti® mostly important 
resonant region. Since ^rr(w) « Crs(‘^) sufficiently 
high frequency w, see Sec. II one can conclude that 
would vanish in the resonant region as well. Note that 
if the latter limit is not yet reached, then one observes a 
contribution at resonance, which corresponds to the error 
^Rg(a;) —This explains our numerical results for 
the realistic systems studied, namely the absence of the 
signal for the ionic liquid and the small but finite contri¬ 
bution to the resonant region observed for the HOD in 
water case, see insets in panels b) and c). Interestingly, 
this line of reasoning suggests that having the frequency 
high can only make the situation worse, as then the error 
would become smaller and the spectral density even more 
underestimated. In fact, considering HOD in D 2 O, where 
a frequency separation between the OH stretch and the 
bath modes was provided, fully confirmed this conjecture, 
did not improve the result (data not shown). Still, the 
low frequency region, lu ^ uj, comes out correctly as it 
follows from this discussion and has been demonstrated 
numerically above. 

Strictly speaking, the obtained spectral density is gen¬ 
erally valid only for the particular (high frequency) mode 
studied, thus making this low frequency region irrelevant 
and the whole rigid bond approach inappropriate. An 
important consequence of the generally dramatic under¬ 
estimation of the spectral density in the resonant region 
is the correspondingly overestimated energy relaxation 
time given in the framework of the Landau-Teller the¬ 
ory nni simply as = 1 /^lp(w). 

One might ask at this point, why the rigid bond 
method had success for the A 2 in A system, which we 
elaborate on in the next section. 


Unexpected success of rigid bond method for A2 in 

A 

Let us consider what would happen if the system stud¬ 
ied was indeed of the CL form. Then Crb(^) = Ccl(^) 
without any approximations. If further the model system 
potential, Vs(a:), was harmonic then ^cl(0 Clp(^) 
would coincide up to a trivial offset. Since Crb( 0 = 
^Lp(i) again up to a constant offset, ^rs( 0 needs not to 
be considered. It suggests that the success of the rigid 
bond method for A 2 in A is based on the possibility that 
this system is indeed of the CL form. 

In fact, Berne et al. [33] have tested this conjecture by 
considering the temperature dependence of the static fric¬ 
tion, = 0) nnd the frequency renormalisation term, 

^Ch{t = 0). They concluded that the temperature depen¬ 
dence observed excludes the possibility of such a direct 
correspondence. However, the temperature was probed 
starting just above the melting point, which might have 
led to a significant change of the properties. Addition- 




1000 2000 3000 4000 3300 3600 3900 
tf)(cm'^) 0) (cm'^) 


FIG. 8: Real parts of the spectral densities for harmonic 
(red) and anharmonic (blue) system potentials are shown in 
the left panels for a) A 2 in A b) HOD in bulk water. Right 
panels zoom on the respective resonant regions. 


ally, we have recently observed that the failure of such 
an indirect check does not necessarily lead to visible dis¬ 
crepancies in the observables, e.g. linear vibrational spec¬ 
tra. [H] Particularly, testing the linearities of the cou¬ 
pling on both the system and bath sides for the very 
same realistic systems suggested that the latter is much 
more important than the former. Performing the same 
test for A 2 in A led to a strongly non-linear coupling on 
the system side and an almost ideally linear one on the 
bath side (data not shown). Therefore, the important 
linearity on the bath side suggests that the system can 
be represented by the CL model. 

Another possible check is provided by the indepen¬ 
dence of the spectral density in the CL model from the 
system potential, see Eq. (0.8). This can be tested by 
computing the memory kernel, ^cl(^)) for both harmonic 
and anharmonic system potentials. It turned out that the 
spectral densities due to (an) harmonic system potentials 
coincide for A 2 in A, see panel a) in Fig. A reference 
calculation performed for HOD in water revealed a strik¬ 
ing discrepancy in the resonant region, panel b) therein. 

These tests considered together indicate that A 2 in A 
can indeed be represented by the CL model. This also 
suggests that testing the dependence of the spectral den¬ 
sity on the system potential is plausible. 


CONCLUSIONS AND OUTLOOK 

In this paper we considered the problem of parametriz¬ 
ing spectral densities within the linear GLE framework. 
We have shown that from the general perspective the 
problem naturally poses itself in the frequency domain, 
and both time and Laplace domains are not intrinsic 
to the problem. We have developed a Fourier-based 
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method and compared its performance to the existing 
time-domain techniques, that is Berne and Harp |44| and 
the widely used rigid bond approach. It turned out that 
our method is extremely robust and provides trustworthy 
results for all systems studied. In contrast, Berne and 
Harp approach suffers from numerical problems apart 
from having general disadvantages of the time domain 
formulations. Surprisingly the rigid bond method, which 
is claimed to be accurate in the high frequency limit, 
turns out to be not appropriate in general, unless the 
system studied is of the Caldeira-Leggett form. Impor¬ 
tantly, we have shown that the spectral densities due to 
the rigid bond method are strongly underestimated in 
the resonant region, which makes the common Landau- 
Teller estimation for the relaxation times questionable. 
A parametrizing scheme beyond the linear GLE regime, 
that is based on Eq. (0.14), has also been developed and 
will be published elsewhere. 

The authors gratefully acknowledge financial support 
by the Deutsche Forschungsgemeinschaft (Sfb 652). 
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